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Abstract 

CFD-based aeroelastic computations are performed to investigate the effect of nonlinear aerodynamics on transonic 
LCO (Limit Cycle Oscillations) characteristics of a two-dimensional supercritical wing with the NLR 7301 section. It 
is found that the presentation of the viscous effects, including turbulence modeling, plays an important role on the 
accurate prediction of shock and LCO; and a small initial perturbation appears to produce large amplitude LCO at 
small mean pitch angle and plunge while a large amplitude initial perturbation produces small (or negligible) 
amplitude LCO at larger mean values. Also addressed in the paper is the issues related to multiblock MP1 (Message 
Passing Interface) parallel computation. 


Introduction 

LCO has been a persistent problem on several current 
lighter aircraft and is generally encountered with 
external store configurations. Denegri [1] provided a 
detailed description of the aircraft/store LCO 
phenomenon. Norton [2] gave an excellent overview of 
LCO of fighter aircraft carrying external stores and its 
sensitivity to store carriage configuration and mass 
properties. 

LCO can be characterized as sustained periodic 
oscillations which neither increase nor decrease in 
amplitude over time for a given flight condition. Using 
an s-domain unsteady aerodynamic model of the 
aircraft and stores, Chen, Sarhaddi and Liu [3] have 
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shown that wing/ store LCO can be a post- flutter 
phenomenon whenever the flutter mode contains low 
unstable damping. This type of flutter mode is called a 
“ hump mode'". Since the aircraft structure usually 
contains structural nonlinearity such as friction 
damping, this amplitude-dependent friction damping 
can suppress the growth of amplitude, thus resulting in 
a steady state oscillation. This is known as the 
nonlinear structural damping (NSD) model of the 
wing/store LCO. Although not thoroughly proven 
through tests or numerical simulations, results of the 
NSD show excellent correlation with flight test LCO 
data of F-16 throughout subsonic and transonic Mach 
numbers. On the other hand, other researchers, notably 
Cunningham and Meijer [4], believe that the wing/store 
LCO is due largely to the transonic shock oscillation 
and shock induced flow separation, called Transonic 
Shock/ Separation (TSS) model. Edwards has suggested 
the TSS model and viscous effects are two major 
factors that cause transonic LCO for wings. He also 
has studied the shock buffet phenomenon in addition to 
transonic LCO [5], It should be noted that, however, 
there is no conflict in the NSD model and the TSS 
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model in that both physical effects may contribute to 
LCO. 

Recent renewed interest in LCO is perhaps motivated 
by the need to further understand the physics of LCO 
and the current advent of CFD methodology in 
aeroelasticity. There are two potential computational 
methods for LCO prediction/investigation: the CFL3D 
code (version 6) [6-8] developed and supported by 
NASA/Langley and the POD/ROM EigenMode 
approach [9] originated by Dowell and Hall of Duke 
University. The former is a conventional time-domain 
CFD method whereas the latter a frequency-domain 
CFD method, using aerodynamic eigenmodes. 

The present study plans to use a CFD time-marching 
method, CFL3D v6, to numerically investigate 
transonic LCO of a supercritical airfoil under a 
plunging/pitching spring-mounting system [10-12], It is 
reasonable to start from investigating a two- 
dimensional LCO case in order to better understand the 
physics of LCO. However, because of the complexity 
of a two-dimensional LCO experimental test, there is 
few experimental data available for comparison. 
According to our best knowledge, the experimental 
work performed by Schewe et. al. [10-12] is perhaps 
the only two-dimensional LCO experimental test 
available in documents. Those test data were 
immediately used by Platzer et. al. to validate their 
thin-layer Navier-Stokes aeroelastic solver [13-14]. 
While the emphasis of [13-14] was on the predictive 
capability of the thin-layer Navier-Stokes aeroelastic 
solver, our emphasis here is to investigate the effect of 
nonlinear aerodynamics on transonic LCO of the 
supercritical airfoil, such as the impact of the viscous 
terms, different turbulence modeling, and the initial 
disturbances. Also addressed in the paper is the issues 
related to multiblock MPI parallel aeroelastic 
computation. 

Numerical Methodology 

The computer code used in this study is CFL3D v6, 
which solves the three-dimensional thin-layer Reynolds 
averaged Navier-Stokes equations with an upwind finite 
volume formulation [6], A two-dimensional problem 
can be calculated by using two identical grid planes, 
created by duplicating the two-dimensional grid. 

The code uses formally third-order upwind-biased 
spatial differencing for the inviscid terms with flux 
limiting in the presence of shocks. Either flux- 
difference splitting or flux-vector splitting is available. 
The flux-difference splitting method of Roe [15] is 
employed in the present computations to obtain fluxes 
at cell faces. Viscous terms are discretized with second- 


order central differencing. There are two types of time 
discretization available in the code. The first-order 
backward time differencing is used for steady 
calculation while the second-order backward time 
differencing with i-TS subiterations is used for static 
and dynamic aeroelastic calculation. Furthermore, grid 
sequencing for steady state and multigrid and local 
pseudo-time stepping for time marching solutions are 
employed. Also available in the code are many 
turbulence models, although here only the Spalart- 
Allmaras model [16] and Baldwin-Lomax model [17] 
with the Degani-Schiff modification have been used. A 
detailed description of the methodology of the code can 
be found in [6], 

One of the important features of the CFL3D code is its 
capability of solving multiple zone grids with one-to- 
one connectivity. Spatial accuracy is maintained at zone 
boundaries, although subiterative updating of boundary 
information is required. Coarse-grained parallelization 
using the MPI protocol can be utilized in multiblock 
computations by solving one or more blocks per 
processor. When there are more blocks than processors, 
optimal performance is achieved by allocating an equal 
number of blocks to each processor. As a result, the 
time required for a CFD-based aeroelastic computation 
can be dramatically reduced. In this paper, both single 
and multiblock MPI parallel aeroelastic computations 
near the onset of flutter LCO are compared with 
experiment and with other computations. Figure 1 
shows a C-type grid with 273 x 93 mesh points around 
the NLR 7301 airfoil that has been divided into eight 
69x47 blocks. This and a single block version of this 
grid are used in the computations to follow. 



Figure 1 Multiblock C-type grid around NLR 7301 
airfoil (eight 69x47 blocks) 


Accordingly the mesh deformation scheme in [7] is 
modified to fit multiblock grids. In [7], the mesh 
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deformation uses a modified spring analogy with solid 
body translation/rotation of the fluid mesh near solid 
surfaces. Initialization of the grid deformation at each 
step is performed using a TFI (Transfmite 
Interpolation) step. The mesh interior is then smoothed 
and grid orientation near boundaries is preserved using 
the modified spring analogy. In the present 

implementation, the subgrid based TFI scheme of [18] 
has been employed for initialization at each time step. 
That scheme uses subgrids consisting of “slave 
vertices” to move both block boundaries and interiors. 
In some instances, in order to achieve an optimal 
division of grid points, it is necessary to place flow 
field block boundaries near a moving solid surface. An 
example of this is shown in Figure 1. The multiblock 
boundary and interior movement scheme allows the 
user to place block boundaries near surfaces as 
necessary for optimal parallelization. Boundaries 
interior to the fluid domain near a surface respond to 
the local surface motion. As the airfoil moves, block 
boundaries move to maintain integrity of block 
interfaces and the airfoil surface. User controlled input 
makes it possible to update the mesh using this 
subgrid/TFI-based scheme alone or to update with an 
initialization using this scheme plus additional 
smoothing steps. These added smoothing steps, the 
number of which can be defined by the user, employ 
the modified spring analogy scheme [7]. In the current 
implementation the spring analogy scheme updating the 
mesh interior is now written in delta formulation so that 
the relative orientation of the original grid is retained. 
The solid body rotation/translation of the fluid grid is 
also now performed near both solid surface and block 
fluid boundaries. 

The time-marching simulation of the aeroelastic 
responses is obtained using the state transition matrix 
solution from t to t+At of the state variable 
representation of the decoupled modal equations [19- 
20], The state transition matrix based scheme is optimal 
in the sense that it is derived from an exact solution of 
the free response of the modal equations. The actual 
scheme uses predictor/corrector steps. The predictor 
step marches the structure using the solution of the 
modal equations at the step n to get the surface 
deflection at the time step n+1. This provides the 
surface shape for a recomputation of the fluid mesh and 
the fluid domain solution at n+1. After a solution of the 
fluid domain involving multiple subiterations, the 
corrector step then solves the modal equations at the 
time step n+1 using the averaged generalized forces at 
n and n+1. 

Because the CFD and CSM meshes usually do not 
match at the interface, CFD/CSM coupling requires a 
surface spline interpolation between the two domains. 


The interpolation of CSM mode shapes to CFD surface 
grid points is done as a preprocessing step. Modal 
deflections at all CFD surface grids are first generated. 
Modal data at these points are then segmented based on 
the splitting of the flow field blocks. Mode shape 
displacements located at CFD surface grid points of 
each segment are used in the integration of the 
generalized modal forces and in the computation of the 
deflection of the deformed surface. The final surface 
deformation at each time step is a linear superposition 
of all the modal deflections. 



Figure 2 Two-degree-of-freedom dynamic model 

The following is an account of our theoretical modeling 
of Schewe’s experiment on transonic flutter of a two- 
dimensional supercritical wing with an NLR7301 airfoil 
section [10-12]. Figure 2 depicts a simplified model of 
the two-degree-of-freedom test set-up. The two- 
dimensional wing has a chord length of 0.3 m (c = 0.3 
m ) and a span of 1 m ( b = 1 m). The pitching spring and 
the plunging spring are attached to the same c/4 
position. The corresponding two-degree-of-freedom 
equation of motion of the set-up reads 

I m h - A l h l l Dh 0 l it I 

L ■ a 4m _ 1 <• j L 0 A _ 1 - ] 

, I K h 0 l h i . i L(t) i 

L 0 K _ ] „ j - 1 M{t) J 

where m h is the total mass ( m h = 26.64 kg), I cl4 is 
the mass moment of inertia about c/4 ( I c/4 = 0.086 kg- 
m), s\ is the static unbalance ( y = 0.378 kg-m), D h 
and D are the damping factors of the plunging motion 
(/?) and the pitching motion (i) respectively (D h = 
82.9 kg/s and D — 0.197 kg-m 2 /(rad-s)), K h and K : 
are the stiffness of the plunging spring and the pitching 
spring respectively ( K h = 1.21 x 10 6 N/m and K : = 
6.68 x 10 3 N-m/rad), and L(i) and M{t) are the 
aerodynamic lift and moment respectively in Newtons. 
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The aeroelastic equations and the CFD grid are 
maintained in dimensional form. To perform the time- 
marching CFD computation in CFL3D v6.0, it is 
necessary to convert Eq. (1) into modal coordinates, 
i.e.: 


! \ 1*11 q , ( 2 ) 

I L J 

where q is the modal coordinate and « is the modal 
matrix of the undamped structure. For this numerical 
example we have 


I - 0.1735 0.1004 
4 " [ 0.9277 3.403 


(3) 


Substituting Eq. (2) into Eq. (1) and pre-multip lying the 
resulting equation by u ' yields 


|/| 


I 2 o h L h 

L 0 



1<M l 


* h 


0 

2 \q- 


- V 


T I L(t) I 
]M(t) j 


(4) 


where t h and 1 t are the undamped natural 
frequencies of the plunging and pitching motions 
respectively ( o h = 205.4 rad/s and o t = 299.5 rad/s), 
l h and l u are the plunging and pitching damping 
ratios respectively (<. h = 0.00648 and l t = 0.00474). 
Note that the off-diagonal terms in the damping matrix 
are assumed to be zero for simplicity. 


Results and Discussions 


The simulated case here is the measurement No. 77 
documented in [12]. As mentioned before, the 
experimental model was a two-dimensional 
supercritical wing with NLR7301 section. The chord 
length of the wing was 0.3 m and the angle of attack 
was 1.28 . The experimental conditions were the free- 
stream Mach number of 0.768 and the Reynolds 
number of 1.727/ 10 6 based on the chord length. A 
transonic LCO in two-degrees-of-ffeedom was found at 
the dynamic pressure of 0.126 bar. The corresponding 
free-stream velocity was 254.7 m/s. The total pressure 
was 0.45 bar. 


The criterion used in [13-14] was to match the 
computed to the measured time-averaged surface 
pressure distribution. 

Time-Averaged Surface Pressure Distribution 

First, an Euler computation is performed on a C-type 
grid with 293x61 points. The best agreement with the 
experimental data is at M= 0.734 and t - - 0.25 . 
Flowever, even for this corrected Mach number and the 
corrected angle of attack, the predicted shock strength 
is stronger than the experimental result and the location 
of the shock is behind the measured one, as shown in 
Figure 3. We have searched all Mach numbers and 
angles of attack. It seems impossible to match both the 
strength and location of the shock with Euler 
computation. Then two viscous computations with 
Baldwin-Lomax turbulence model and Spalart- 
Allmaras turbulence model are performed on a C-type 
grid with 293x93 points. The corrected Mach number is 
found to be 0.748 for both models while the corrected 
angle of attack is -0.02 for the Baldwin-Lomax model 
and 0.15 for the Spalart-Allmaras model. Figure 3 
indicates that both viscous results have a closer 
agreement with the experimental data, especially for 
shock strength and location, clearly showing that 
viscous effects are important for the accurate prediction 
of the shock. 



Figure 3 Time-averaged surface pressure 
distribution (Aerodynamics only) 


As shown in [13-14], because of the relatively large 
chord length of the airfoil with respect to the wind 
tunnel test section (1 m x 1 m), both the freestream 
Mach number and the angle of attack need to be 
corrected to take into account wind tunnel wall effects. 


The above computed results are obtained from steady 
CFD analyses of a rigid two-dimensional wing surface. 
With taking into account the effects of elastic 
responses, the corrected angle of attack becomes -0.1 
for Euler computation, 0.15 for the Baldwin-Lomax 
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model, and 0.32 for the Spalart-Allmaras model. The 
corrected Mach number does not change for both Euler 
computation and the Baldwin-Lomax model but shifts 
to 0.75 for the Spalart-Allmaras model. These 
parameters are used in the following dynamic 
aeroelastic computations. Again, as shown in Figure 4, 
the result from Euler computation overpredicts the 
shock strength. The match between the computed 
pressures of the viscous solutions and experiment is 
slightly less near the shock. 



Figure 4 Time-averaged surface pressure 
distribution (static aeroelastic computation) 


Performance and Time Step Convergence 

The static aeroelastic case is used to compare run times 
between the single and 8 block/MPI parallel 
computations. Starting from a steady state a time- 
accurate aeroelastic solution is marched for 800 time 
steps at 5 multigrid subiterations each. The Spalart- 
Allmaras turbulence model is used in each case. 

Table 1 gives performance data for the various grids, 
run modes, and code configurations. The first two rows 
in the table represent sequential computation of the 1 
and 8 block versions of the grid. The 8-block grid takes 
28% less CPU time than the 1 -block grid, apparently 
due to better caching of the smaller sized blocks. The 
parallel computation with 8 processors runs 7.8 times 
faster than the same run sequentially. The last two rows 
present the computational effort required for the spring 
analogy smoothing step and the combined TFI and 
aeroelasticity. Three iterations of the spring analogy 
scheme add about 1.5% computing time to the solution. 
The TFI grid movement and the aeroelastic 
computation add about 10% to the run time. 


Table 1. Cost/P erformance 


Grid/comp. 

type 

TFI 

Spr. Anal. 
Smoothing 

CPU 

Time/ 

Processor 

(seconds) 

1 block/seq 

yes 

yes 

8800 

8 block/seq 

yes 

yes 

5430 

8 block/MPI 

yes 

yes 

690 

8 block/MPI 

yes 

no 

680 

8 block/MPI 1 

no 

no 

610 


+ also without aeroelasticity 


Table 2 shows solution behavior for the single block 
grid for successive time step sizes. These are 
computed at M„=0.753, a=0.6 and a dynamic pressure 
of 0.126 bar, using the Spalart-Allmaras model. 
Initialization is accomplished with a static aeroelastic 
solution, followed by an initial perturbation of the 
dynamic simulation of -.00114 in the velocity of the 
second mode. The time step sizes give 80, 250, and 800 
time steps per pitch/plunge cycle. Nine subiterations per 
time step are used. Columns two and three are half 
amplitudes of the fully developed LCO plunge and 
pitch, while columns four and five arc the plunge/pitch 
frequencies. Frequencies and amplitudes are computed 
based on a sampling of the last 8-10 cycles of motion; 
from this sampling the data appears to be nearly 
converged at the smallest time step. At the largest time 
step, even after 100 cycles, the amplitude slowly 
continued to grow while at the two smaller time steps 
the amplitude fully converged to LCO. In all of the 
remaining dynamic computations, the smallest time 
step size with nine multigrid subiterations is used. 


Table 2, Time step convergence, 1 block grid 


At 

a h (m) 

a«(Deg.) 

C0h(Hz) 

Wa(Hz) 

0.128 

.0112 

3.72 

34.5 

34.5 

0.040 

.00917 

3.24 

34.3 

34.3 

0.0125 

.00899 

3.17 

34.3 

34.3 


The speed increase offered by computing in parallel is 
appealing. There are trade offs of course when using 
coarse gram parallelization. Depending on the block 
splitting and problem, the multiblock computations can 
require substantially more subiterations. This fact is 
most evident in the problem at hand. Figures 5 and 6 
show a comparison between single block and 
multiblock limit cycle plunge and pitch. As an aside, it 
must be stated that after several hundred cycles the 
amplitude and frequency of the multiblock pitch and 
plunge had settled out and converged to values virtually 
identical to that of the single block grid (a h = 0.00887 
m, a „= 3.13 , 0 )|, = 0 )„ = 34.3 Hz). This would not be the 


5 


American Institute of Aeronautics and Astronautics 





case if the multiblock aeroelastic coupling and 
integration were not consistent with that of the single 
block configuration. Yet as shown in figures 5 and 6, 
the early evolution of the multiblock LCO is much 
different than that of the single block. With these results 
casting some doubt on the convergence of the 
multiblock CFD for this LCO problem, the final 
computations are completed with the single block grid. 


0.02 r 


0.01 h 


h (m)o 


- 0.01 


>.02 





1 - 5 1 1 1 ! i 


8 block 
1 block 

I 1 1 1 1 I 1 I 1 1 I 1 ; I 1 I 


5 1.75 2 2.25 2.5 

t(s) 


Figure 5 Plunge vs. time (N-S, small perturbation, 
Mo = 0.753, a = 0.6°) 
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Figure 6 Pitch vs. time (N-S, small perturbation, 
Mo = 0.753, a = 0.6°) 


Effect of Viscosity and Turbulence Model 

Figures 7 and 8 suggest that without the viscous terms, 
the predicted plunge and pitch modes appear never 
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towards LCO. The effect of viscosity is to limit the 
amplitude of the flutter LCO. This is made clear by 
figures 5 and 6 and by figures 9 and 10. Figures 9 and 
10 present dynamic plunge/pitch responses to a small 
perturbation (-.00114 in the velocity of the second 
mode) using the Baldwin-Lomax and the Spalart- 
Allmaras models. At the time of writing, the Spalart- 
Allmaras simulation of Figures 9-10 had not reached a 
limit cycle whereas the Baldwin-Lomax reached large 
amplitude limit cycle quite rapidly. Plunge and pitch 
amplitudes of the B-L result are around 0.004 m and 2 
respectively. It is clear that the turbulence model 
significantly alters the nature of the solution. The S-A 
results in each of the cases shown in figures 9-10 and 
figures 5-6 appear to be approaching a fixed point LCO 
while the B-L results of figures 9-10 are LCO but 
chaotic in nature. 



Figure 7 Plunge vs. time (Euler, small perturbation, 
Mo = 0.734, a = -0.1°) 



Figure 8 Pitch vs. time (Euler, small perturbation. 
Mo = 0.734, a = -0.1°) 
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Figure 9 Plunge vs. time (N-S, small perturbation, 
M n = 0.750, □=0.32^ (S-A), 

Mo = 0.748, □= 0.15 s (B-L)) 


(Deg.)0 



response to the small perturbations both yielded mean 
plunge and pitch very near their static equilibrium 
positions. In response to the large perturbation, both the 
amplitude and mean plunge and pitch had offset 
significantly. Mean plunge had offset to values of 
0.004 m (B-L) and 0.003 m (S-A). Mean pitch angles 
had offset to values of 3 degrees (B-L) and 1.9 degrees 
(S-A). Amplitudes of both turbulence model 
simulations are much smaller in response to the large 
disturbance and both show evidence of continuing 
chaos. The amplitudes of plunge are around 0.002 m 
(B-L) and 0.001 m (S-A). Amplitudes of pitch are 
around 1 degree (B-L) and 0.2 degree (S-A). 




Figure 11 Plunge and pitch vs. time (B-L, large 
perturbation, Mu = 0.748, □=0.15°) 



Figure 10 Pitch vs. time (N-S, small perturbation, 
M n = 0.750, □= 0.32° (S-A), 

Mu = 0.748, □= 0.15° (B-L)) 

Effect of Perturbation Size 

The effect of initial perturbation size is studied by 
repeating earlier simulations with both turbulence 
models, but with a large initial perturbation (-.114 in 
the velocity of the second mode). Results are shown in 
figures 11 and 12. After the expected very large 
transients passed, the solutions using the two turbulence 
models appear to have reached LCO solutions. In 



Figure 12 Plunge and pitch vs. time (S-A, large 
perturbation, M D = 0.750, □= 0.32°) 
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The plunge and pitch amplitudes computed here with 
the single and multiblock grids, due to small 
perturbation, are generally in the range of those 
computed elsewhere. In particular, the presently 
computed pitch amplitude due to the S-A model (a D = 
3.17 degrees) agrees very well with other results, while 
the plunge amplitude appears to be much lower (a h = 

0.009 m computed here). [13] These pitch amplitudes 
are much larger than the pitch amplitudes of the 
experiment. (,2-.4 degrees) The mean pitch angle at .7 
degree computed here compares with a mean angle of 
1.28 degrees in the experiment. 

The amplitudes of responses due to the large 
perturbation actually are much closer to those of 
experiment. The pitch amplitude for the S-A model 
(around .2 degrees) and mean pitch (around 2 degrees) 
agrees quite well with experiment. It must be cautioned 
that the data under discussion (i.e. of figure 12) may not 
have reached LCO, and may in fact be very slowly 
diminishing in amplitude. What can be said is that the 
small initial perturbation produced large amplitude 
LCO at small mean pitch angle and plunge while a 
large amplitude initial perturbation produced small (or 
negligible) amplitudes at larger mean values. Future 
computations with intermediate values of perturbation 
may produce somewhat larger amplitude (than the 
small or negligible) and lower mean pitch (than the 
large mean values) that produce a better match with 
experiment. 

A final comment about the B-L results must be made. 
The chaotic nature of the B-L results must be taken 
with a certain amount of caution. To what extent the 
chaotic nature of the LCO is due to numerical as 
opposed to physically based reasons is not certain. 
What is remarkable is that both turbulence models 
produce a similar reduction in LCO amplitudes as a 
result of the large initial disturbance. 

Conclusions 

Viscous effects are found very important for the 
prediction of both shock and LCO. Special attention 
should be paid to the numerical representation of 
viscous effects, including turbulence modeling. A small 
initial perturbation appears to produce large amplitude 
LCO at small mean pitch angle and plunge while a 
large amplitude initial perturbation produces small (or 
negligible) amplitude LCO at larger mean values. 
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